home *** CD-ROM | disk | FTP | other *** search
/ Night Owl 6 / Night Owl's Shareware - PDSI-006 - Night Owl Corp (1990).iso / 039a / d3d.zip / CABLE.C < prev    next >
Text File  |  1991-08-30  |  7KB  |  257 lines
  1. /* CABLE: This program reads a file that represents a space curve, and
  2.    writes a file that represents a cable.  The circular cable section
  3.    will be a regular polygon with n vertices; it approximates a circle
  4.    with radius R.  Both n and R are read from the keyboard.  The program
  5.    is to be linked together with the module TRAFO, in which the functions
  6.    'initrotate ' and 'rotate' are defined.
  7. */
  8.  
  9. #include <stdio.h>
  10. #include <math.h>
  11. #include <alloc.h>
  12. #include <process.h>
  13.  
  14. void  initrotate( double a1, double a2, double a3,
  15.                   double v1, double v2, double v3, double alpha);
  16. void rotate(double x, double y, double z,
  17.       double *px1, double *py1, double *pz1);
  18.  
  19. int zero(double x);
  20.  
  21. double *getdouble(int n);
  22.  
  23. double *enlarge(double *px, int N);
  24.  
  25. void ermes(char *s);
  26.  
  27. main()
  28. {
  29.    char infil[30], outfil[30];
  30.     int i, n, j, m, jn, jn0, k, tablesize;
  31.    double R, *x, *y, *z, xC0, yC0, zC0, xC1, yC1, zC1, xC2, yC2, zC2,
  32.       *xC, *yC, *zC, a, b, c, d, rx, ry, rz, pi, theta, Len, *getdouble(),
  33.       *enlarge(), xA, xB, yA, yB, zA, zB, dx, dy, dz, d0, c1, c2, c0,
  34.       xM, yM, zM, e1, e2, e0, denom, lambda, mu, xP, yP, zP, xAP, yAP, zAP,
  35.       xBP, yBP, zBP, v1, v2, v3, cosphi, phi;
  36.    FILE *fpin, *fpout;
  37.    printf("Input file: ");
  38.     scanf("%s", infil);
  39.       if((fpin = fopen(infil, "r")) == NULL)
  40.       ermes("File not available");
  41.    if
  42.     (fscanf(fpin, "%*d %lf %lf %lf", &xC0, &yC0, &zC0) != 3 ||
  43.      fscanf(fpin, "%*d %lf %lf %lf", &xC1, &yC1, &zC1) != 3 ||
  44.      fscanf(fpin, "%*d %lf %lf %lf", &xC2, &yC2, &zC2) != 3)
  45.       ermes("Input file incorrect");
  46.  
  47.    printf("Output file: "); scanf("%s", outfil);
  48.    if((fpout = fopen(outfil, "w")) == NULL) 
  49.       ermes("Problem with output file");
  50.  
  51.    printf("How many points on each circle? "); scanf("%d", &n);
  52.     printf("Radius: "); scanf("%lf", &R);
  53.  
  54.    a=xC2 - xC0; 
  55.    b=yC2 - yC0; 
  56.    c=zC2 - zC0; 
  57.    d = a*xC1 + b*yC1 + c*zC1;
  58.  
  59.    /* First circle has center (xC1, yC1, zC1), radius R, and lies
  60.       in the plane ax + by + cz = d  */
  61.    x = getdouble(n);
  62.    y = getdouble(n);
  63.    z = getdouble(n);
  64.  
  65.    if(zero(a) && zero(b)) {
  66.       rx =  0;
  67.       ry =  c;
  68.       rz = -b;
  69.    }
  70.    else {
  71.       rx =  b;
  72.       ry = -a;
  73.       rz =  0;
  74.    }
  75.  
  76.    Len = sqrt(rx*rx + ry*ry + rz*rz);
  77.    rx /= Len;
  78.    ry /= Len;
  79.    rz /= Len;
  80.    /* (rx, ry, rz) is a unit vector perpendicular to (a, b, c) */
  81.    x[0] = xC1 + rx*R;
  82.    y[0] = yC1 + ry*R;
  83.    z[0] = zC1 + rz*R;
  84.    pi = 4.0 * atan(1.0);
  85.    theta = 2.0*pi/n;
  86.  
  87. /* Computation of n points on the first circle */
  88.    initrotate(xC1, yC1, zC1, a, b, c, theta);
  89.    for(i=1; i<n; i++)
  90.       rotate(x[i-1], y[i-1], z[i-1], x+i, y+i, z+i); /* x+i == &x[i] */
  91.  
  92. /* Count the number of circles (number of points -1 from input file) */
  93.    m = 2;
  94.    while(fscanf(fpin, "%*d %lf %lf %lf", &xC0, &yC0, &zC0) == 3)
  95.       m++;
  96.    tablesize = m*n;
  97.    x = enlarge(x, tablesize);
  98.    y = enlarge(y, tablesize);
  99.    z = enlarge(z, tablesize);
  100.    rewind(fpin);
  101.    fscanf(fpin, "%*d %lf %lf %lf", &xC0, &yC0, &zC0); /* Skip 1 */
  102.  
  103.    xC = getdouble(m);      /* get pointers to memory tracts */
  104.    yC = getdouble(m);
  105.    zC = getdouble(m);
  106.  
  107.    for(j=0; j<m; j++)
  108.       fscanf(fpin, "%*d %lf %lf %lf", xC+j, yC+j, zC+j);
  109.    fclose(fpin);
  110.  
  111.    /* (xC[0], yC[0], zC[0]) is now the center of the given circle,
  112.       lying in plane ax + by + cz = d, and with radius R. The n
  113.       relevant points on this circle have already been computed;
  114.       their coordinates are
  115.             x[0], y[0], z[0], ..., x[n-1],y[n-1],z[n-1].
  116.       The other m-1 circles will be derived from this first one by means
  117.       of rotations.
  118.    */
  119.  
  120.     for(j=1; j<m; j++) {
  121.         jn  = j*n;
  122.         jn0 = jn - n;
  123.  
  124.       xA = xC[j-1];
  125.       yA = yC[j-1];
  126.       zA = zC[j-1];
  127.  
  128.       xB = xC[j];
  129.       yB = yC[j];
  130.       zB = zC[j];
  131.  
  132.       dx = xB - xA;
  133.       dy = yB - yA;
  134.       dz = zB - zA;
  135.  
  136.       c1 = a*a  + b*b  + c*c;
  137.         c2 = a*dx + b*dy + c*dz;
  138.         c0 = d - a*xA - b*yA - c*zA;
  139.  
  140.       xM = 0.5*(xA + xB);
  141.       yM = 0.5*(yA + yB);
  142.       zM = 0.5*(zA + zB);
  143.  
  144.       d0 = dx*xM + dy*yM + dz*zM;
  145.  
  146.       e1 = dx*a  + dy*b  + dz*c;
  147.       e2 = dx*dx + dy*dy + dz*dz;
  148.       e0 = d0 - dx*xA - dy*yA - dz*zA;
  149.  
  150.       denom = c1*e2 - c2*e1;
  151.  
  152.       if(fabs(denom) < 1.e-12) {
  153.          /* Direction does not change.
  154.             Instead of using a point P indefinately far away,
  155.             we perform a simple translation: */
  156.  
  157.          for(i=0; i<n; i++) {
  158.             x[jn+i] = x[jn0+i] +dx;
  159.             y[jn+i] = y[jn0+i] +dy;
  160.             z[jn+i] = z[jn0+i] +dz;
  161.          }
  162.       }
  163.       else {
  164.          /* Direction changes.
  165.             The polygon will be rotated through the angle phi
  166.             about the vector v passing through the point P. */
  167.  
  168.             lambda = (c0*e2 - c2*e0)/denom;
  169.             mu = (c1*e0 - c0*e1)/denom;
  170.          xP = xA + lambda*a + mu*dx;
  171.          yP = yA + lambda*b + mu*dy;
  172.          zP = zA + lambda*c + mu*dz;
  173.  
  174.       /* Point P(of intersection of three planes) is center of rotation */
  175.             
  176.          xAP = xA - xP;
  177.          yAP = yA - yP;
  178.          zAP = zA - zP;
  179.  
  180.          xBP = xB - xP;
  181.          yBP = yB - yP;
  182.          zBP = zB - zP;
  183.  
  184.          v1 = yAP*zBP - yBP*zAP;
  185.          v2 = zAP*xBP - zBP*xAP;
  186.          v3 = xAP*yBP - xBP*yAP;
  187.  
  188.          /* (v1, v2, v3) is the direction of the axis of rotation */
  189.  
  190.          cosphi = (xAP*xBP + yAP*yBP + zAP*zBP)/
  191.                   sqrt((xAP*xAP + yAP*yAP + zAP*zAP)*
  192.                        (xBP*xBP + yBP*yBP + zBP*zBP));
  193.          phi = (cosphi == 0 ? 0.5*pi : 
  194.                atan(sqrt(1.0 - cosphi*cosphi)/cosphi));
  195.          /* phi is the angle of rotation */
  196.          initrotate(xP, yP, zP, v1, v2, v3, phi);
  197.          for(i=0; i<n; i++)
  198.             rotate(x[jn0+i], y[jn0+i], z[jn0+i], x+jn+i, y+jn+i, z+jn+i); 
  199.          initrotate(0.0, 0.0, 0.0, v1, v2, v3, phi);
  200.          rotate(a, b, c, &a, &b, &c);
  201.       }
  202.       d = a*xC[j] + b*yC[j] + c*zC[j];
  203.    }
  204.    for(k=0; k<m*n; k++)
  205.         fprintf(fpout, "%d %f %f %f\n", k+1, x[k], y[k], z[k]);
  206.  
  207.    fprintf(fpout, "Faces:\n");
  208.  
  209.    for(k=n-1; k>=0; k--)
  210.       fprintf(fpout, " %d", k+1);
  211.    fprintf(fpout, ".\n");
  212.  
  213.     for(k=(m-1)*n; k<m*n; k++)
  214.       fprintf(fpout, " %d", k+1);
  215.    fprintf(fpout, ".\n");
  216.  
  217.    for(j=1; j<m; j++) {
  218.       jn = j*n + 1;
  219.       jn0 = jn - n;
  220.       for(i=0; i<n-1; i++) {
  221.          fprintf(fpout, " %d %d %d.\n", jn+i+1, -(jn0+i), jn0+i+1);
  222.          fprintf(fpout, " %d %d %d.\n", jn0+i, -(jn+i+1), jn+i);
  223.       }
  224.       fprintf(fpout, " %d %d %d.\n", jn, -(jn0+n-1), jn0);
  225.       fprintf(fpout, " %d %d %d.\n", jn0+n-1, -jn, jn+n-1);
  226.    }
  227.    fclose(fpout);
  228. }
  229.  
  230. int zero(double x)
  231. {
  232.    return fabs(x) < 1.e-5;
  233. }
  234.  
  235. double *getdouble(int n)
  236. {
  237.    double *p;
  238.     if((p = (double *)malloc(n*sizeof(double))) == NULL)
  239.       ermes("Not enough memory");
  240.    return p;
  241. }
  242.  
  243. double *enlarge(double *px, int N)
  244. {
  245.     if((px = (double *)realloc(px, N*sizeof(double))) == NULL)
  246.       ermes("Not enough memory");
  247.    return px;
  248. }
  249.  
  250. void ermes(char *s)
  251. {
  252.    printf(s);
  253.    exit(1);
  254. }
  255.  
  256. 
  257.